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Abstract 

We will present results of a numerical integration of a maximally sliced Schwarzschild black hole using 
a smooth lattice method. The results show no signs of any instability forming during the evolutions to 
t — lOOOrn. The principle features of our method are i) the use of a lattice to record the geometry, ii) the 
use of local Riemann normal coordinates to apply the 1+1 ADM equations to the lattice and iii) the use of 
the Bianchi identities to assist in the computation of the curvatures. No other special techniques are used. 
The evolution is unconstrained and the ADM equations are used in their standard form. 

PACS numbers: 04.25.Dm, 04.60.Nc, 02.70.-c 



1. Introduction 
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J ■ Recent studies [1-7] have shown that the stabihty of numerical integrations of the Einstein 

^' field equations can depend on the formulation of the evolution equations. Subtle changes in 

the structure of the evolution equations have been shown to have a dramatic effect on the long 

/\ • term stability of the integrations. These are relatively new investigations and thus at present 

d • there is no precise mathematical explanation as to what is the root cause of the instabilities 

or how best they can be avoided or minimised. What we have at present is a growing set 

of examples which suggests that the standard ADM evolution equations may not be the 

most suitable equations for numerical relativity. Consequently many people are looking at 

alternative formulations such as the hyperbolic formulations of Einstein's equations [8-10] 

and the conformal ADM equations of Shibita and Nakamura [6] and Baumgarte and Shapiro 

[3]. 

One alternative is the smooth lattice approach which we presented in two earlier papers 
[11,12]. This is a method which uses a lattice similar to that used in the Regge calculus 
but differing significantly in the way the field equations are imposed on the lattice. In the 
smooth lattice method we employ a series of local Riemann normal coordinates in which the 
connection vanishes at the origin of each such frame. Collectively these frames enable us to 
obtain point estimates of the curvatures in terms of the lattice data (in particular the leg 
lengths). The upshot is that the 3+1 ADM equations can be applied directly to the lattice. 



This is clearly a radically different approach to that normally used in numerical relativity. It 
is thus interesting to explore its stability properties against those for traditional techniques. 

In our first paper [11] we showed how the smooth lattice method could be used to obtain the 
initial data for a Schwarzschild spacetime. In the second paper [12] we showed how the 3+1 
ADM equations could be applied to a lattice using the Kasner spacetime as a test case. In 
both papers the results were very encouraging. In this paper we return to the Schwarzschild 
spacetime, this time to study the stability of its evolution in a maximal guage. 

Studies of a maximally sliced Schwarzschild spacetimes in spherical symmetry were quite 
popular some years ago (see for example [13,14]). These studies showed that the evolutions 
were invariably unstable. The source of the instability was attributed to the stretching of the 
grid as grid points nearer the black hole were drawn into the black hole quicker than those 
further out. With the consequent loss of resolution the estimates for the derivatives were 
seriously in error and the non-linear feedback in the equations quickly drove the solution 
into exponential overfiow. 

We will repeat these calculations using our smooth lattice method so that we can address 
the simple questions : When will the loss of resolution become apparent, what impact will 
it have on the subsequent evolution and will it trigger an unstable evolution? 

We should point out that this paper is not an attempt to revive the use of maximal slicing as 
a preferred slicing condition. We are, instead, using it solely as a test of the smooth lattice 
method. 

We will try to stay as close as possible to the earlier work of Bernstein, Hobill and Smarr 
[14]. Thus we shall not be using any of the modern techniques, such as apparent horizon 
boundary conditions [15], conformal differentiation [16] and using the Hamiltonian constraint 
to stabilise the maximal slicing equation [13]. We will, however, employ some techniques of 
our own, in particular we will use a lattice to record the metric, we will cover the lattice with 
a series of local frames and we will use the Bianchi identities in computing the curvatures. 

As this method may be unfamiliar to many readers we have included more details of the 
derivations than might normally be included. However, to spare the reader we have relegated 
the bulk of the derivations to a (large) appendix. The sections that precede the appendix 
contain all of the important results with few derivations. We begin with a brief review of 
the smooth lattice method, followed by a description of the particular lattice used for the 
Schwarzschild spacetime. We then present the details of the 1+1 ADM equations, the results 
of the integrations and finally we review what we have found. 
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2. Riemann Normal Coordinates and Smooth Lattices 

The concept of a smooth lattice was first introduced by Brewin [11,12] as a method by which 
a discrete lattice could be coupled to a family of Riemann normal frames in such a way as 
to provide smooth estimates for the metric and curvature on the lattice. 

The principle features of the smooth lattice method are 

o The lattice is a finite collection of vertices and legs connected in any suitable fashion 
(e.g. a simplicial lattice). 

o The data recorded on the lattice is purely geometrical, such as the leg lengths and angles 
between pairs of legs. 

o To each vertex their is assigned a small neighbourhood in which local Riemann normal 
coordinates are employed. Each such neighbourhood is called a computational cell. For 
a simplicial lattice the computational cell can be chosen as the set of simplicies attached 
to the central vertex. 

o The legs of a computational cell are taken as geodesic segments of a locally smooth 
metric. 

o The leg lengths are small compared to the curvature length scales, that is, RL^ « 1 
where R and L typical values for the curvature and leg lengths respectively. 

There are a number of equivalent definitions of Riemann normal coordinates [17,18,19,20]. 
One definition has them to be the coordinates which, for a given point P, the geodesies 
through P are all of the form x^{s) = sa^ where s is an afiine parameter and the a^ are 
constants (i.e. they are "straight" lines). The coordinates of a point Q near P are then taken 
as x^{s) for the geodesic that joins P to Q. Clearly this restricts the region in which these 
coordinates can be used to that for which there is a unique geodesic which joins P to Q. In 
larger regions it is possible that pairs of geodesies may cross and thus the coordinates at the 
intersection would not be unique. This problem is avoided in a smooth lattice by requiring 
the leg lengths to be small when compared to the curvature length scales, i.e. RL? « 1. 

An equivalent definition of Riemann normal coordinates is that the connection and its sym- 
metric first derivative vanish at P, 

o = r^«/3 (2.1) 

= T''aP,X + r^^A,a + r^Aa,/? (2.2) 

In a Riemann normal frame it is relatively easy to show that the metric, when expanded 
about the point P, takes the form 

g^^A^) = g^^. - \R^.aupx'^xP + 0(x3) (2.3) 
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where g^u are constants that can always be chosen to be diag(l, 1, 1). This prescription does 
not uniquely determine the coordinates as rotations around P are still allowed. This guage 
freedom can be used to orient the axes to preferred edges of the lattice. 

From this form of the metric it is relatively easy to compute various geometrical quantities 
such as the geodesic length Lij between two vertices i and j 

1% = g^.Ax'l^Ax'^j - ^R^a.px'lx^x'^x'^ + ^(e^) (2.4) 

and the angle 6i subtended at the vertex i in the geodesic triangle (ijk) 

2UjL,k cos 9, = 4 + Ll - L% - ii?^„,^A4 A^. Aa^ffcAxf, + 0{e'>) (2.5) 

In these equations x'^ are the coordinates of vertex i, Axf- = xf — x'^ and e is a typical 
(small) length scale. See [11] for a derivation of these equations. 

The principle advantages of using a smooth lattice method over other lattice methods are 

o The metric is smooth and differentiable. 

o Point estimates of the curvatures are easy to compute. 

o The equivalence principle is explicitly used in each computational cell. 

o The vanishing of the connection at a vertex greatly simplifies many equations (i.e. co- 
variant differentiation reduces to partial differentiation). 

3. The lattice 

The smooth lattice method could equally be applied to the 4-dimensional spacetime or, in 
a 3+1 ADM context, to each spacelike Cauchy surface. We will adopt the second approach 
simply because it reduces the complexity of the bookkeeping (a 3-d computational cell has 
far less complexity than its 4-d counterpart). The 3+1 ADM equations will be used to evolve 
the 3-d smooth lattice. 

For a Schwarzschild spacetime it is reasonable to choose the Cauchy surface to be a lattice 
built upon concentric 2-spheres. Each two sphere could be subdivided into a set of cells (e.g. 
triangles) using the same pattern on every 2-sphere. The successive 2-spheres can be joined 
by short radial legs connecting pairs of similar vertices. This construction is not unique as 
it does allow for a creeping rotation to occur between successive 2-spheres. This can be 
eliminated by demanding that the radial legs be perpendicular to the 2-spheres. 

As part of the smooth lattice approach we require each leg in the lattice to be a geodesic 
segment of the 3-metric. Thus, as the radial legs are required to be normal to each 2-sphere, 
we see that each sequence of connected radial legs also forms a global geodesic of the 3-metric. 
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What flexibility do we have in choosing the leg lengths? Since the space must by spherically 
symmetric we can only adjust the overall scale of each two sphere and the distance between 
successive 2-spheres. These two pieces of information can be recorded by specifying the 
leg lengths between a pair of 2-spheres, L^z, and the typical length within a 2-sphere, Lxx- 
By this means we can reduce the complexity of a full 3-dimensional lattice to a simple 
2-dimensional ladder as indicated in flgure 2. 

3.1. The computational cell 

Each Riemann normal frame was chosen to cover the region between three consecutive 2- 
spheres. Thus each Riemann normal frame will include the three rungs L^^i ^%xi ^xx ^^^ 
the two struts L^^ and L~^. We also chose to align each Riemann normal frame so that i) 
the ladder was confined to the a;2;-plane, ii) the z-axis coincided with the familiar radial axis 
and iii) the 2;-axis threaded the mid-points of each rung (see figures (1,2,3)). It follows that 
the coordinates of the vertices must be of the following form 

Table 1. Coordinates for the vertices in Figure 3 



Vertex {x, y, z) 


Vertex (a;, y, z) 


Vertex (x, y, z) 


r ( u, 0, V-) 


1 ( M°, 0, 0) 

2 {-u\ 0, 0) 


1" ( n\ 0, v^) 
2^ {-u\ 0, v^) 



with f < and v" > 0. Note that the origin of the Riemann normal frame has been located 
(by a translation in the radial direction) so that the z coordinate of vertices 1 and 2 are zero. 

We also require that successive pairs of struts form a global radial geodesic, that is, there 
can be no kink at the vertex where the radial struts meet. Thus we demand that 



IT 



d^ + d~ (3.1.1) 



which we call the geodesic constraint. From these five leg lengths and one geodesic constraint 
we need to compute the curvatures and all of the coordinates (see section 8.1 for details). 

3.2. The Riemann curvatures 

For a spherically symmetric space there are just two algebraically independent curvature 
terms, Rxyxy and Rxzxz = Ryzyz (the equality arises from the rotational symmetry around 
the z-axis). In Appendix 8.1 we show that the system of equations for the leg lengths and 
the geodesic constraint can be reduced to a single equation for Rxzxz 

9 / f^ — r° r~ — r° \ 

n / XX ^xx I ^xx ^xx iiD ro (o n ^\ 

U = j+ . J- T^ \ TZ j + Rxzxz^xx (3.2.1) 



zz ' zz \ zz zz 
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which in turn is seen to be a finite difference approximation (on a non-uniform lattice) to 
the geodesic deviation equation for nearby radial geodesies, namely 

(PL 

U = , 2 I" -^xzxz-l^xx [O.Z.I) 

The remaining smooth lattice equations serve only to determine the coordinates u°, -u"*", etc. 
Clearly we need one more equation in order to compute the second curvature Rxyxy- This 
could be obtained by introducing extra structure into the lattice, such as the diagonal braces 
joining pairs of vertices on a 2-sphere. In fact this structure already exists - when we first 
spoke of the lattice we imagined each 2-sphere to be fully triangulated. Only later latter did 
we choose one of the legs on which to build our ladder. We could simply choose the collection 
of triangles attached to a particular vertex as a base on which to build a more sophisticated 
lattice with one ladder built over each leg of each triangle. This lattice would contain two 
classes of ladders - those sharing the common radial geodesic (that generated by the central 
vertex) and a chain of ladders forming a cylinder around the common radial geodesic. This 
lattice would therefore contain two classes of rungs - one for each class of ladder. In principle 
this lattice should allow us to compute both curvatures, Rxzxz and Rxyxy However, there 
is a further complexity in that we do not know, a priori, the relationship between the leg 
lengths of the two classes of rungs. One might be tempted to set the lengths of each class of 
rung in the asymptotically fiat region (i.e. far from the throat) and to then impose the same 
ratios (between the two classes) for all of the rungs back down to the throat. This would 
be correct if the rungs of the ladders where geodesic segments of each 2-sphere. However, 
the rungs are geodesic segments of the full 3-dimensional metric and thus their ratios will 
change with distance from the black hole's throat. Rather than pursue a solution to this 
problem we chose instead to retain our simple (one ladder) lattice and to employ the Bianchi 
identities to compute the second curvature Rxyxy For our spherically symmetric space we 
can show (see Appendix 8.7) that there is only one non-trivial Bianchi identity 

r. _ [Lxx^xyxy] - [Lxx^xyxy] 1 / + , po \ (-^xx) ~ (-^xx) Ci. O 'i.\ 

^~ T+ r) {-^XZXZ ' -'^XZXZJ T+ yO.Z.O) 

^ZZ ^ ^zz 

which is is a simple forward finite difference approximation to the continuum equation 

d [L^^Rxyxy) dL^^ (ooA\ 

U = ; iixzxz — ; — [6.1.^) 

dz dz 

in which z is the proper distance measured along the radial axis from the throat. 

We solve the coupled equations (3.2.1,3.2.3) for the curvatures given all the leg lengths and a 
suitable initial value for Rxyxy on the inner boundary of the lattice (see section 5 for details). 
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4. The 1+1 ADM equations 



4.1. The evolution equations 

In an earlier paper [12] we showed how the 3+1 ADM evolution equations may be applied to 
any lattice. For the present problem with zero shift and drift the evolution equations may 
be written (see Appendix 8.2) as 



dL 



2NK.,Axf,Ax'^. (4.1.1) 



^^ -'-M— y— *j 



d 



- (/f^.Axf.Axr,) = (-ATj^, + N {R^, + KK^, - 2K^aK)) A^A^- (4.1.2) 



which when applied to our lattice (see Appendix 8.2) leads to 

^ = -NK,,L^^ (4.1.3) 

dt 

^ = -NK^^L,, (4.1.4) 

dK 

-^ = -N,^^ + N {R^^ + KK^^) (4.1.5) 

^ = -iV,z. + N {R,, + KK,,) (4.1.6) 

wnere A = ZJ\xx ~r J^zzi J^xx -^ -'t x^x ^ J^xyxy + J^xzxz SHQ ix^z -^ -'t z^z ^ ^J^xzxz- -l-ilG 

partial derivatives of the lapse function N^xx and N^zz can be evaluated using the techniques 
in Appendix 8.5, leading to 

1 dLxxdN /. -, „x 

N^xx = 7 — 1^-r (^•^•^) 

d^N 
N„ = ^ (4.1.8) 



4.2. The constraint equations 

For the Schwarzschild spacetime there are only two non-trivial constraint equations, the 
Hamiltonian constraint 

= R + K^-K>"'K^„ (4.2.1) 

and the momentum constraint 

= K\^- K/\, (4.2.2) 
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where K = K^ ^. These equations are readily adapted to the smooth lattice (see Appendix 
8.6) leading to 

U = -ttxyxy ~r ^Hxzxz ~\~ -t^xx ' ^^xx^zz [^.Z.Oj 

= ^ JLxxKxx) _ ^^^'E^ (42.4) 

dz dz 



5. Numerical methods 

5.1. The initial data 

The initial data for the lattice consists of the leg lengths Lxx-, Lzz and the extrinsic curvatures 
Kxx) Kzz- These can be freely chosen subject to the two constraints (4.2.3) and (4.2.4). For 
a time symmetric slice we must have = Kxx = Kzz and consequently the momentum 
constraint (4.2.4) is identically satisfied. This leaves the Hamiltonian constraint, which now 
takes the simple form, 

U ^ -tlxyxy + Z-ttxZXZ [O.l.lj 

as the one equation to constrain Lxx and Lzz- Following Bernstein, Hobill and Smarr [14] 
we chose to set Lzz while computing the Lxx as a solution of the Hamiltonian constraint. 

The radial leg lengths Lzz where set as follows. A stretched grid of isotropic radial coordinates 
were defined by 

r,- = -e(^^^) (5.1.2) 

with j = on the inner boundary (the throat) and j = N on the outer boundary. The 
parameters A and A^ were chosen by Bernstein, Hobill and Smarr so that the outer boundary 
was at r ^ 200m. They chose A = 6/A^ and A^ = 200 while for our production runs we 
chose A^ = 800. The Lzz were then chosen as 

1+1 / rn\'2 



for j = to j = A^ - 1. 

The Lxx were set by re-arranging the coupled system (3.2.1,3.2.3) and (5.1.1) in the form of 
a radial integration. Starting from the throat and working outwards, 

1+ _ I 

^xx = ^xx + T^r [L^^ — L^^) + -tL^^ [L^-. + L^-,) [LxxRxyxy) (5.1.4) 

P+ = f?° I \^xx) ~ [^xx) \ (^^ ^\ 

-^xyxy -^xyxy \ ^,^2 \+ (t2 \° I yo.i.o) 

\ ^ K-'^xx) " K-^xa 



■'XXI 



■^xzxz cy^xyxy [O.l.Xj) 



At the throat we imposed reflection symmetry by setting L^^. = L~^, L^^ = L~^. We also 
chose L°^ = m/10 (the equations are linear in L°^ and thus this choice is not crucial) and 
Rxyxy = l/(4m^) = —2R°^^^ (which we obtained from the analytic solution). Finally, we set 
m = 1. With this information we applied the above equations to generate all of the lattice 
data from the throat out to the outer boundary. 

Is it reasonable to be using the analytic solution to assist us in setting the initial data? The 
only information borrowed from the analytic solution are the Lzz in each cell and Rxyxy at 
the throat. Where we chose to locate the successive Riemann cells is up to us, that is, we 
are free to choose the L^z as we see fit. This is identical to the freedom to choose the lapse 
function when evolving the initial data. Thus this use of the analytic solution is not crucial 
and is made only to allow us to make direct comparisons with Bernstein, Hobill and Smarr. 
But what of the choice of curvature Rxyxy'^ Notice that the analytic solution depends on 
just one parameter, the mass m. Our lattice initial data also depends on just one parameter, 
the value of Rxyxy at the throat. Thus whatever choice we make for Rxyxy we are in effect 
choosing the ADM mass of our numerical spacetime. We could make some other choice for 
Rxyxy and then later determine the ADM mass for our numerical spacetime. To take this 
approach would be tedious and thus we chose to take the easier option where we set the ADM 
mass at the outset. Note that we make no further use of the analytic solution throughout 
the subsequent evolution. 

5.2. The boundary conditions 

The standard boundary conditions for Schwarzschild initial data are that the inner boundary 
is reflection symmetric and that in the distant regions the data is asymptotically flat. The 
reflection symmetry can be imposed by extending the lattice so as to have a computational 
cell that straddles the throat and then to demand that the two halves of this cell be mirror 
copies of each other. With this extra computational cell we can apply any of the lattice 
equations at the throat. 

For a reflection symmetric throat we demand L^^, = L~^ and L'l^ = L~^ at the throat. 
This condition was imposed throughout the evolution by setting L~^ and L~^ equal to their 
updated counterparts L'^x ^^^ ^tz ^^ each stage of the integration (i.e. within each of the 
four steps of the 4th order Runge-Kutta). 

At the outer boundary we extended the lattice by half a cell so that we could apply the 
evolution equations to the data associated with that cell. The data for the extra half cell 
were obtained by cubic extrapolation from the interior. The only data that needed to be 
extrapolated was Kxx, N, Nzz and Rxzxz for the grid centred scheme (see section (5.4)) and 
Rxzxz for the standard scheme. In both schemes we set = dL^.^/dt. 
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5.3. The lapse function 

A maximally sliced spacetime is defined to be a spacetime for which K = everywhere. 
This is normally imposed by setting i^ = on the initial Cauchy surface and then setting 
dK/dt = throughout the evolution. This condition leads, through the standard ADM 
evolution equations, to the following elliptic equation for the lapse function 

= \/^N-RN (5.3.1) 

Using the results of section 8.5 we can write this as 

= — ^ + - ^ -RN (5.3.2) 

dz^ Lxx dz dz 

which in turn can be applied to the lattice by replacing each of the derivatives by finite 
difference approximations on a non-uniform lattice (see Appendix 8.8) 

The common wisdom is that to obtain a stable evolution, the Hamiltonian constraint should 
be used to eliminate the curvatures terms in the above equation. However, we found (see 
section 6) that the evolution was stable without this modification. 

The boundary conditions of reflection symmetry at the throat, dN/dz = 0, and asymptotic 
flatness limj-^oo A^ = 1 were prescribed on our (finite) lattice as 

inner boundary : = A^ — A^~ (5.3.3) 

outer boundary : 1 = A^° (5.3.4) 

The coupled system (5.3.2-5.3.4) were solved by three iterations of a shooting method. In 
each shot we made a guess A^~ for A^~ on the inner boundary and then used a Thomas 
algorithm to solve a modified system of equations consisting of A^~ = A^~ at the inner 
boundary, the main equation (5.3.2) everywhere on the lattice (except the outer boundary) 
and the outer boundary equation 1 = A^°. The two initial guesses for A^~ were taken as 
A^~ = and A^~ = 1. After each shot, an estimate for dN/dz at the inner boundary was 
formed. 

Since the equation (5.3.2) is linear in A^ it is possible to form a linear combination of the two 
solutions (one from each shot) so as to satisfy both the inner and outer boundary conditions. 
This leads to the third and final guess 

with the subscripts denoting the shot number. This guess was used for the final shot for the 
lapse on the lattice. 
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5.4. Time stepping 

In each Riemann normal frame we can use the evolution equations to compute the time 
derivatives of the lattice data. However some of the lattice data are shared between neigh- 
bouring frames, L^z for example, and thus we will obtain multiple estimates for their time 
derivatives. How then should we compute a single estimate for the time derivatives of the 
lattice data? Though there are many schemes which could be imposed, we chose two related 
schemes. 

In the first scheme, which we will refer to as the standard scheme, we formed simple averages 
over all of the time derivatives for each component of the lattice data. For example, as each 
leg joins two vertices there will be two estimates available for its time derivative, one from 
each vertex. For some legs, such as Lxx-, the two vertices are in equivalent frames (i.e. the 
two frames are at the same distance from the throat) and thus yield identical values for 
the time derivatives. Thus, for these legs it was sufficient to compute just the one time 
derivative. However, for legs such as Lzz-, the two Riemann frames are distinct and both 
time derivatives must be computed. In the standard scheme no averaging was performed for 
the time derivatives of the K^p. 

In the second scheme we treat Lzz and Kzz as if they are defined at the centre of their 
associated legs. All other data, Lxx^ Kxx, Rxyxy, Rxzxz, N, Nxx, Nzz are taken to be defined 
on the vertices. When the time derivatives are evaluated, interpolations of the data are 
required to assemble all terms at the appropriate point. We take the average of the vertex 
data when interpolating from a pair of vertices to the centre of the leg. However, the 
interpolation of the data from the centres of the legs to the vertices is slightly more involved. 
For a function / defined at the centre of the radial legs the interpolation to the common 
vertex of two successive radial legs is given by 

/ = J, ]j- {Llzf + L:J^) (5.4.1) 

^zz "•" ^zz 

This formula was only applied to the Kzz terms when computing the time derivatives for 
Lxx and Kxx- The final step was to take averages as per the standard scheme. This method 
will be referred to as the grid centred scheme. 

One could ask why the Lxx and Kxx where not given a similar centred treatment. The 
answer is that, because of the rotational symmetry of the lattice, the interpolation would 
produce values identical to that obtained by assuming the Lxx were based on the vertices. 
So for simplicity we choose to take Lxx and Kxx as based on the vertices in the grid centred 
scheme. 

For all of our production runs we chose a 4-th order Runge-Kutta scheme with a fixed time 
step of 5t = 0.01. The evolution equations (4.1.3-4.1.6) were treated as a fully coupled 
system of equations. 
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In each of the four steps of a single Runge-Kutta cycle we would first compute all of the 
Rxzxz using (3.2.1). From (3.2.3) we would then compute the Rxyxy, from the throat to the 
outer boundary, using the Hamiltonian constraint (4.2.3) as a boundary condition for Rxyxy 
at the throat. Then followed the computation of the lapse function, the time derivatives 
and their averages before the partial updates were made. This would be repeated for the 
remaining steps in the Runge-Kutta cycle. 

6. Results 

The results of our integrations, using the grid centred scheme on a grid with 800 radial legs, 
are presented in figures (5-13). Notice that all of the curves are smooth and show no signs 
of any instabilities to t = IOOOtti. We have also marked the location of the apparent horizon 
on each curve with a diamond (the procedure for locating the apparent horizon is described 
later in section 6.5). This clearly shows that the horizon propagates smoothly and with 
almost constant speed across the grid. It is also clear that majority of the dynamics occurs 
within a very narrow region straddling the horizon. 

The standard scheme produced curves that were qualitatively similar to those from the grid 
centred scheme. In particular they showed no signs of any numerical instabilities out to 
t = 1000m. They did, however, show non-trivial quantitative differences in the later stages 
of the evolution - at t = 1000m the standard scheme's grid had stretched out to 739m as 
opposed to 826m for the grid centred scheme. The grid centred scheme also showed better 
constancy in the area of the apparent horizon (see section (6.5)), with only a 4% change in 
area from t = to t = 100m as opposed to a 13% change for the standard scheme. In the 
later stages of the evolution we can expect that the accuracy of the solution will be degraded 
due to the loss of resolution near the apparent horizon. This problem will occur with all 
numerical methods that do not provide for special treatment over the highly dynamical 
regions (e.g. an adaptive grid refinement scheme such as that due to Berger and Oliger [21]). 

There are a number of simple checks that have been used by many others [13,14,16] to check 
the quality of the numerical solution. These include monitoring the constraints, the growth 
of the apparent horizon, the convergence of the maximal slices to r = 3m/2 and the so called 
collapse of the lapse. We will discuss each of these checks in turn. 

6.1. Constraints 

Ideally the constraints should remain bounded throughout the evolution. In practice though 
the constraints do drift away from their initial values. The plots in figures (11,12) do show 
an initial growth in the constraints but after t ^ 100m the constraints remain bounded. 
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6.2. The r = 3m/2 limit 

It has been shown by Estabrook etal [22] that the maximal slices of the Schwarzschild 
spacetime converge to the limit surface r = 3m/2 where in this case r is the standard 
Schwarzschild radial coordinate. We can use this as a check on our solution. Those parts 
of the grid for which the lapse has collapsed will be frozen on r = 3m/2. Thus the lattice 
data should be constant across those sections. This can be seen in the early stages of the 
evolution where inside the apparent horizon all of the data are (approximately) constant. 
For the curvature terms we can estimate what these constants should be. Notice that the 
geometry inside the apparent horizon is that of a cylinder, i.e. S'^ x R. The scalar 3-curvature 
will be just that of the 2-curvature of the 2-dimensional cross sections (i.e of a 2-sphere of 
radius r = 3m/2). Thus we must have R = 2/r^. We also know that R = 2(Rxyxy + '2Rxzxz) 
and as L^^ is constant along the cylinder we find from (3.2.2) that Rxzxz = 0. Thus we deduce 
that Rxyxy = 4/(9m^) = OA/m^ which agrees well with our numerical value of 0.44445/?72^ 
during the early evolution, from t = to t = 100m. However, we loose agreement in the 
long term evolution, t = 100m to t = 1000m. This is caused by a loss of resolution near 
the apparent horizon due to grid stretching. To show this, we plotted Rxyxy at t = 100m 
for four different grid resolutions, with 100, 200, 400 and 800 grid points, see figure (14). 
This clearly shows that the grid resolution has a significant impact on the accuracy of the 
solution. 

6.3. Geodesic slicing 

It is well known [23] that for a geodesicly sliced Schwarzschild spacetime (i.e. setting the 
lapse equal to 1) the throat will collide with the r = singularity at coordinate time nm. 
This provides a simple test - run the code and see when it crashes (its also a curious way 
to compute tt!). We ran the code for 100, 200, 400 and 800 grid points and found that the 
code crashed within two time steps of t = 3.14 (with a time step of 0.01). 

6.4. Collapse of the lapse 

Beig [24] has shown that for a maximally sliced Schwarzschild spacetime, the lapse at the 
throat will die exponentially with time, 

iV ~ /5e~"* + C (e-2"*) ast^oo (6.4.1) 



where 



a = —^ ^ 0.54433 (6.4.2) 

3v6 

(3 = -^ exp ( ^ ) ^ 0.83725 (6.4.3) 



3^2 VsVe. 

^ = ^ ln(54v^ - 72) - 2 In [ '^^'^ \ ~ _o.21815 (6.4.4) 

4 ' I 9^6 -22 ' 
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See also the earlier works by Estabrook and others [22, 25, 26, 27]. This gives us another test 
where a plot of In A^ versus t should be a straight line. This we have done in figures (15) and 
for t = 10m to t = 100m we get a very straight line indeed. Using a standard least squares 
method we found a = 0.5474 and f3 = 1.011. Though our estimate of a agrees well with 
Beig's result, our estimate for f3 is not so good. This is probably due to a number of factors 
such as the use of of an inexact outer boundary condition, the problems of grid stretching and 
numerical error. The numerical error in fitting a straight line could be significant. Notice 
that for < t < 100m the vertical intercept in figure (15) is very small (approximately 
0.01) relative to the vertical range (approximately 60). Thus any small errors, either in the 
data or in fitting the line may produce large relative errors in the vertical intercept, namely 
ln/3. To demonstrate this point we chose to re-compute a and f3 by constraining the curve 
to pass through the data point at t = 10m. Thus we applied the least squares method to 
IniV — In A^(IO) = a{t — 10) and found a = 0.5429 and f3 = 0.8112. This is an improvement 
over our earlier estimates. 

At later times, t = 100m to t = 1000m the line does bend slightly. This is not surprising 
since we know that for these times we have lost accuracy due to grid stretching. Note that 
at t = 1000m we have N ~ 10~ at the throat - the evolution really has been halted. 

6.5. Apparent horizon 

An apparent horizon is defined as a closed 2-surface with zero divergence of its outward 
pointing null vectors [28]. For the Schwarzschild spacetime the apparent horizon must be 
a 2-sphere whose area must be constant when Lie dragged in the direction of the outward 
pointing null vectors. Thus 

= CuA=^ (6.5.1) 

au 

where A is the area of the 2-sphere and d/du is the outward pointing null vector. For our 
lattice we may put 

(6.5.2) 

(6.5.3) 

where /c is a pure constant, z is the proper distance measured along the radial axis and A^ is 
the lapse function. Thus from = dA/du we obtain the apparent horizon equation 

= ^ - L^xK^:, (6.5.4) 

dz 

Our aim is to solve this equation for two things i) the location of the apparent horizon and 
ii) the value of L'^^ at that location. We do so as follows. First we convert the derivative in 
Lxx into a finite difference approximation using the methods of section 8.8. Then we scan 



-14- 



A-- 


- ^^xx 


d 
du 


1 d d 

~ Ndt^ dz 



the grid while monitoring the value of the right hand side of (6.5.4) which we denote by 
the function Q{z). When we find two successive grid points for which Q changes sign, we 
stop the scan and then use linear interpolation to predict the true location of the apparent 
horizon (i.e. z such that = Q) and any data at that point (e.g. L^^ and A^). 

Since there is no gravitational radiation in the Schwarzschild spacetime we know that the area 
of the apparent horizon must remain constant throughout the evolution. We have plotted 
L^^ = A/k in figure (16). For the case of 800 grid points we see that the area varies by less 
than 4% for the first 100m of the evolution. However at later times, around t = 1000m, the 
error has blown out by 700%. This is very large indeed. We do not believe that this is an 
error in our code but rather a very severe consequence of the loss of resolution due to grid 
stretching. 

Anninos etal. [16] found that for their 1-dimensional code, with a constant Ar = 0.1m and 
an outer boundary at r = 130m, that at t = 25m their apparent horizon mass (defined by 
Mah = {A/16n) ' ) had grown by about 4% (see their figures 15 and 16). When we repeated 
their calculations, on their grid but using our equations, we found an error of 3% for the 
standard scheme and 1% for the grid centred scheme. In a related work, Anninos etal. [15] 
found that for a grid with 400 points the error in the apparent horizon mass was about 25% 
at t = 100m. For our code, with 400 grid points, we found errors of 12% for the standard 
scheme and 6% for the grid centred scheme. 

We should mention here that setting A = kL^^ is not strictly correct because the rungs of the 
ladder are not geodesic segments of the 2-sphere but rather of the surrounding 3-dimensional 
space. However if the L^^ are small relative to the area of the apparent horizon then there 
will be little error in using these legs as approximate geodesies on the apparent horizon. 

7. Discussion 

The results we have presented are very encouraging. Our most important observation is that 
the evolution is very stable out to t = 1000m. In contrast, the best results using traditional 
methods (15,16) developed fatal instabilities by about t = 100. Thus we are forced to ask 
the obvious question : What is it about our method that, in this instance, gave us a stable 
evolution? The features that distinguish our method from traditional methods is that we 
use a lattice to record the metric, we cover the lattice with a series of Riemann normal 
frames and we use the Bianchi identities to assist in the computation of the curvatures. It 
is premature at this stage to identify which, if any, of these ingredients is crucial to the 
stability, but we can review each of their roles in the method. 

Riemann normal coordinates. These were chosen not simply because they allowed us 
to extract curvatures from the lattice but rather for their relationship with the Einstein 
equivalence principle. In its full 4-dimensional setting the equivalence principle is equiva- 
lent to having the 4-connection vanish in the freely falling frame. Using such frames as a 
computational tool must surely bring some advantage to the computations. However, we 
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are actually using the Riemann normal coordinates only for the 3-geometry and thus this 
argument is not so strong. On the other hand, a zero connection does greatly simplify many 
of the computations, for example, covariant differentiation reduces to simple partial differ- 
entiation. Furthermore, as the source terms in the ADM equations have a different structure 
to that found in Bernstein, Hobill and Smarr (as a typical example), we can expect that the 
stability properties may differ from those for traditional methods. 

The lattice. This is essential for it provides the structure by which the local Riemann normal 
frames can be connected together to form a global coordinate atlas for the spacetime. Each 
pair of adjacent local Riemann normal frames has a non trivial overlap and through the 
scalar data they share on the lattice the transformation form one frame to another is well 
defined. 

The Bianchi identities. We do not consider the Bianchi identities as central to our method 
for they were introduced only to overcome a limitation in our highly simplified lattice. We 
could test their role by redoing our calculations on a more sophisticated lattice. But we 
should point out that recent work by Christodoulou and Klainerman [29] place great impor- 
tance on the Bianchi identities in their proof of long term stability for weak initial data. 

Another option that may explain the stability is that the method may lack sufficient accuracy, 
due to dispersion or truncation errors, that the sharp peaks needed to trigger an unstable 
mode are never resolved. However, our results were always at least as good as those obtained 
by others (15, 16) which tends to discount this option. 

Clearly more work is required by applying this method to other more challenging spacetimes. 
We shall report on these calculations soon. 

8. Appendix 

8.1. The Riemann normal frame 

It is a simple matter to substitute the coordinates listed in table 1 into the smooth lattice 
equations (2.4) and the geodesic constraint (3.1.1). This leads directly to the following 
equations 

(L-,)' = {2u-f - ^R,,^, {2u-v-f (8.1.2) 

{l:,)' = {2u^f - li?,,,, (2«V)' (8.1.3) 

[L;,]' = {v-f + {u- -u°f- \r,,,, {u%-f (8.1.4) 

{Ll,f = {v'f + {u' -u°f- ii?,,,, (..V)' (8.1.5) 
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Q = ^(u°-u^- ii?,.,,, {u°v^f\ +^(u°-u-- ii?,,,, {n°v-f\ (8.1.6) 

For a given choice of Rxzxz the first five equations can be solved for the coordinates u° , u^ 
etc. The last equation is the geodesic constraint in the form = L°^(cos^"'" + cos6'"). This 
equation will constrain the choice of the curvature Rxzxz- 

Though these equations could be solved using a Newton- Raphson method it would be better 
if we could find an explicit solution. The simple trick to achieving this hinges on the fact that 
these equations serve only as an approximation to the true continuum metric and curvatures 
and are valid only when the domain of the Riemann normal frame is small compared with 
the curvature lengths scales. Thus it is sufficient to solve the leg length equations by a 
perturbation expansion around fiat space. 

Consider the first five equations and let Rxzxz = 0{^) with e taken as our expansion param- 
eter. Then the leading order solution is simply 

(2t.°)o = L°, (8.1.7) 

(2«")o = ^- (8-1-8) 

(2«')o = ^'x (8.1.9) 

{v-)l={L-zzf-{n-uX (8.1.10) 

The next level of approximation is obtained by substituting these back into the above equa- 
tions and solving once again for the five coordinates. The result is 

{2u°\=Ll, (8.1.12) 

_ 1 _ _ 2 

(2m )^= L^^ + -RxzxzL^^ {v )q (8.1.13) 

(2«^)i = ^L + \r.zxzLI, [vX (8.1.14) 

{v-)\ = {L:^ - {u- -u°)l- \Rxzxz {u°v-)l (8.1.15) 

{v^^ = {Ll,f - K -n°)\- \Rxzxz {u°vX (8.1.16) 



This process could be continued to obtain higher order approximations but that would be a 
waste of effort as the smooth lattice equations are valid only to linear terms in the curvatures, 
that is to linear terms in e. 

Substituting these approximations into the geodesic constraint and retaining just the terms 
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linear in e leads to 



r+ _|_ T- 

zz '^ zz 



T+ _ TO 

T + 

^zz 



+ 



L- 



-^xx 



~r rixzxz^T 



xzxz-^xx 



XX 



T + 

^zz 



L: 



-R. 




(8.1.17) 



This is easily seen to be the finite difference approximation to the differential equation 










' -tixzxz-^xx 



J-'xxJ^xzxz 



dz 



l.l.li 



which we recognise, apart from the last term, to be the standard geodesic deviation equation 
applied to the two radial geodesies. By making the simple change of scale Lxx -^ ^L^x we 
can see that for A « 1 (i.e. for very short leg lengths) the last term is insignificant compared 
to the remaining terms (at a fixed position on the radial axis). This is to be expected since 
for the two radial geodesies to be nearly parallel everywhere we must have dLxx/dz « 1. 

Upon deleting this term from the discrete equations we obtain 







^zz 



T + 

XX 



TO 

XX 



l: 



TO 

XX 



^ zz \ ^zz ^zz 



+ tixzxzi-'xx 



(3.2.1) 



as our basic equation from which we can compute the curvature Rxzxz- 

Incidentally this requirement that the two radial geodesies be nearly parallel everywhere also 

_ 2 — 2 — — 

shows that {L^^) >> [u — u°)q and thus to leading order [v )o = ~ {Lzz)- Proceeding in 
this fashion we find the following estimates for the coordinates. 



[2u )-^ — L^^ 

_ 1 _ _ 2 

(2u ) ^ = L^^ + -RxzxzLxx [Lzz) 

\ )l ~ XX + "^Rxzxz^xx y^zz) 
\V )^= 1/^2 — —RxzxzLzz [Lxx) 
[V )^ = L^^ — —RxzxzLzz [Lxx) 



(8.1.19) 
(8.1.20) 

(8.1.21) 

(8.1.22) 

(8.1.23) 



These estimates will be of use later in Appendix 8.2. 
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8.2. The 1+1 ADM evolution equations 

In the paper by Brewin [12] the 3+1 ADM evolution equations for a lattice where given as 



rft2 



2 (|(iV^M-)) Axf^.A^ + Qy (8.2.1^ 



where Qij represented the terms involving the shift and drift vectors. We shall immediately 
set Qij = since this corresponds to our lattice which has zero shift and drift. 

For many reasons (in particular, for ease of numerical integration) it is customary to express 
the evolution equations as a system of first order equations. This will be the main focus of 
this section. 

In [12] two important coordinate frames were used. The first was the Riemann normal frame 
for a given computational cell in one Cauchy surface Sq. These coordinates were denoted 
by x^. Some of these coordinates could be freely chosen (by aligning the coordinate axes) 
on some of the vertices of the cell (in particular, at the origin of the cell) . For the remaining 
vertices the coordinates must be computed during the solution of the lattice equations 

-'-'ij Q^J-v^-^ij^-^ij '^J^^.avfi^i^i^j^ j [o.Z.Z) 

Since the leg lengths are expected to be functions of time we must also expect that in the 
Riemann normal frame the vertex coordinates x^ will also be functions of time. Thus in this 
frame we may choose the shift and drift to be zero at the origin of the cell but must accept 
non-zero values elsewhere. 

The second frame used in [12] was not a Riemann normal frame but one in which the shift 
and drift vectors where set to zero everywhere. These coordinates were denoted by x"^. On 
the initial Cauchy surface the two coordinate frames are identical, x^ = x"^ on Sq, but for 
future times the coordinates will differ. This establishes a time dependent transformation 
between the two coordinate frames. 

Consider for the moment the shadow frame with coordinates x"^. For the case of zero shift 
and drift the standard definition of the extrinsic curvature is 



dt "^'"^'^ 



-2NK'' (8.2.3) 



The leg lengths can be estimated from 

4 = a^M^^^lj (8-2-4) 

But in this frame = dx^^/dt and thus we have 



y 



-2NK"Ax''^Ax'!^ (8.2.5) 



dt '"' '^ '^ 
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which when evaluated on Sq, where the two frames coincide, we have 

dLl 

-^ = -2NK^,Axf^Ax'^j (4.1.1) 

Now consider the time derivative of this last pair of equations. These may be written as 

d^Ll .2^^^^//|.^^//. _ 2^K'Ax';^Ax1^ (8.2.6) 



y 



dfi dt y '^ dt ^"^ '^ '^ 

and 



dt^ 



-2N- (ir^^A^A^.) - 2—K,,Ax>t^Ax\^ (8.2.7) 



which if we now combine with the standard ADM equation 



""" -- -N\^v + N {R^, + KK^, - 2K^^K^) (8.2.8) 



dt 

leads directly to 

^ (ir^,Aa;f^.A<^.) = (-iV,^, + N {R^^ + KK^, - 2Kf,aK)) AxfjAx'^^ (4.1.2) 

on Sq. We have written A^,^jy as opposed to N\ni, since in the Riemann normal frame 
the connection vanishes at the origin and thus A^,^iy = ^lujy on the central vertex of the 
computational cell. 

This completes our stated aim - to derive a pair of first order evolution equations - (4.1.1) and 
(4.1.2). These will now be applied to the lattice. To do so first requires equations (2.4,2.5) 
to be solved for the coordinates a;^ for each vertex in the computational cell. This was done 
in section 8.1 and lead to equations (8.1.19-8.1.23). Substituting these in (4.1.1,4.1.2) and 
retaining only the leading order terms we find 

^ = -NK^^L^^ (4.1.3) 

"^^^ ---NK,,L,, (AAA) 



dt 



-N^x + N{R^^ + KK^^) (4.1.5) 



dt 

dK 

—^ = -N,,, + N {R,, + KK,,) (4.1.6) 

dt 
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8.3. The ADM evolution equations : An alternative derivation 

There is another way in which all of the above equations can be easily developed. The idea 
is to begin with the usual ansatz for a spherically symmetric space such as 

rfs2 _ -N'^{r, t)dt^ + A^ir, t)dr^ + B{r, tf{de^ + sin^ ed(f)'^) (8.3.1) 

and to build a lattice in the 6 = 7r/2 plane by assigning leg lengths according to 

L^a: = 5A0 (8.3.2) 

L,^ = AAr (8.3.3) 

The functions A and B would be evaluated at the centre of each leg while Ar and A(j) would 
be chosen as some suitably small numbers. 

The standard ADM equations, for a metric with zero shift, are 



dgjii, 

i 
dK, 



2NK^^ (8.3.4) 



dt 
^" - -Ni., + N {R^^ + KK^^ - 2K^aK''^) (8.3.5) 



We now ask the simple question: What form do these equations take when applied to a 
lattice? The answer is also simple. Since there is no shift vector, the coordinates of each 
vertex remain constant throughout the evolution, thus = dAxfJdt and we therefore have 

^ ((?^,Axf^.A4) = -2NK^,Ax'ljAx'^j (8.3.6) 

^ (ir^^Axf^. A^) = (-iV|^^ + N {R^, + KK^^ - 2K^aK'',)) A<^. A^ (8.3.7) 

These differ from the smooth lattice equations (4.1.3-4.1.6) only in how the terms N\i^^, and 
R^avp are computed. It is easy to show that for the above metric 

N\rr = Nrr ~ jA^rNr (8.3.8) 

iV|^^ = N^^^ + ^B^rN^r (8.3.9) 

Rrcl>r(j} = -J {A^rB^r — AB^rr) (8.3.10) 

Re^e<t>={^ {A^-{B,rf) (8.3.11) 

= 2B,rRr6r6 - BA^ (^^\ (8.3.12) 
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B' .. 



These equations can be adapted to the lattice by making a coordinate transformation into 
a local Riemann normal frame 



8.3.13) 
8.3.14) 

8.3.15) 



R 



dz 


= Adr 


dx 


= Bd(p 


Nzz-- 


j^2 \rr 


N,ccx = 


- 52^I0<^ 


^xzxz ^ 


1 o 


J^xyxy - 


- ^^Re^e(^ 



8.3.16) 
8.3.17) 
8.3.18) 



and then eliminating, where possible, A and B in favour of Lxx and L^z- This leads directly 
to the smooth lattice equations (3.2.1, 3.2.3, 4.1.7, 4.1.8). The exception is the equation 
for Ro(f)0(f, for which there is no direct counterpart in the smooth lattice equations. For this 
quantity we find 

1 (^^2_(^dL^^ I ^gg^^g) 



R 



xyxy 



h^„ 



dz 



This equation can not be used on this lattice for two reasons. First, there is no clear method 
for determining the parameter A0 from the lattice data L^x and L^z- Second, this equation 
simply does not arise as a consequence of the basic smooth lattice equations (2.4, 2.5). This 
problem could be overcome by employing a different lattice. For example, if we chose a 
lattice in which each 2-sphere was fully triangulated then we could reasonably expect that 
both curvatures could be computed from the smooth lattice equations without reference to 
the Bianchi identities. On such a lattice we should also be able to compute A0. 

This alternative derivation is useful not only in giving us confidence that we have the correct 
lattice equations but it also gives us a technique for quickly adapting a continuum equation 
directly to the lattice. Indeed we could well have chosen this as the primary method by 
which to develop our equations. 

8.4. Covariant differentiation 

Since the connection vanishes in Riemann Normal Coordinates we have, for any vector field, 
that 

^% = <^ (8.4.1) 

at the origin of the RNC cell. 

There are probably many ways in which the partial derivatives could be evaluated, however, 
in this section we shall focus on a method based on coordinate transformations and finite 
differences. 
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The basis of our approach is to import values of the vector field from neighbouring RNC 
cells, by simple coordinate transformations such as rotations and translations, and to then use 
this data in a finite difference approximation. This is exactly the same as parallel transport 
since, once again, the connection vanishes, and thus to leading order we can construct the 
transformations as if we were in flat space. 

We will demonstrate this approach for a spherically symmetric 3-geometry. Consider a 
typical RNC cell with its axes oriented as per figure 3. This cell will be denoted by a° 
while four of its six immediate neighbours will be denoted by ct"'", a^ , being the left and right 
neighbours of c", and o""'", a~ being the cells above and below a° along the radial axis. The 
remaining to cells lie along the y-axis relative to o"° and we will have no need of these cells 
in the following calculations. We will use these superscripts to denote quantities which are 
defined relative to the cell with the same superscript. 

Suppose we have a spherically symmetric vector field v^. Then each cell on a 2-sphere will 
share exactly the same values for v'^ in their respective local RNC frames, that is 



,/^i 



,Mr 



ao 



(8A.2) 



for each i,j. We will now import these values into the central cell a". Consider first the 
cells a° and a^. The RNC coordinates of these frames are related by a rotation (to align 
the directions of their axes) and a translation (to align their origins) which we will write as 



X' 



^^ = A^'^x^° + B^ 



(8.4.3) 



The situation is depicted in figure 4. 

The components A'^i, can be assembled into a rotation matrix 



[M 



cosA6' sinA6' 

1 

- sin A9 cos AO 



i + Ae 






r 








-1 






(8.4.4) 



where A^ is a (small) angle of rotation and where n is taken as the row index. The approx- 
imation that A^ is small is imposed so that the we have an accurate approximation to the 
continuum geometry. By inspection of figure 3 we see that A^ = ALxx/Az. 

Using the standard tensor transformation laws we find that the values of v'^''^ in the RNC 
frame of a° to be 



/fir 



Thus for cells o""*" and o'-'^ we find 



dx'^ 



Af'^v' 



(8.4.5) 



./XT 



v'y" = V 






(8.4.6) 

(8.4.7) 



-23- 



v'^"" = v'° 


- Mv^° 


^1x1 ^ ^xo 


- Aev'° 


v'y^ = vy° 




^izl ^ ^zo 


+ Aev''° 



(8.4.8) 

(8.4.9) 

(8.4.10) 

(8.4.11) 

Since these two cells reside on opposite sides of the yz— plane we can use the above to form 
a centred finite difference approximation for v^^x at the origin of a°, namely 

v^'x= . (8.4.12) 

'Ax ^ ^ 

where Ax is, to leading order, the distance between the origins of cells a''^ and a"*". This is 
easily seen to be 2Lxx, leading to 

v^;.=v^,.= -L^^- (8.4.13) 

vy.x = vy,x = o (8.4.14) 

v'.^ =v'x = --^^^t;"° (8.4.15) 

'"" '"" Lxx Az ^ ' 

The same idea can be applied to derivatives in the y direction, leading to 

(8.4.16) 

(8.4.17) 

(8.4.18) 

which clearly could also have been derived by simple symmetry arguments. The z— derivatives 
are very easy to calculate. Since the RNC frames of the three cells (T°, o"~ and o"""" are related 
by simple translations along the z-axis, the coordinate transformations are trivial and lead 
directly to 

Au^ 
v^.^=v^^ = ^ (8.4.19) 

'^ Az ^ ' 



V'^.y = V\y = 







" ■,y — " ,y — 




1 ALxx ^^zo 
J-'xx ^Z 


v'-y = V\y = 


— 


1 ALxx^_yo 

T A ^ 



\z — ^ ,z 



Ayy 



(8.4.20) 



v':z = ^^. = ^ (8.4.21) 

It is useful to make one small change to the above. We will replace the finite difference 
approximations by their continuum limits (i.e. A \-^ d) to simplify the presentations in the 
following sections. 
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8.5. The Laplacian 

Let be a function which is constant on each 2-sphere of a spherically symmetric space. If 
we put v^ = 0'^ then at the origin of each RNC cell we have v^ = 0^^ since Qf^i, = diag(l, 1, 1 
and = r^Q,^ at the origin of each RNC frame. Since is constant on each 2-sphere, we 
must have, 

<P,x = 0, 0,j, = O, <^,^ = ^ (8-5-1 

on the z— axis. Using this and the results of the previous section we find that, on the z— axis 

(8.5.2 



',xy — V ,y — 
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\xz = V"" ,z = -3- = 

dz 

dvV 
, - .;y - _ - 



^,xx — V ^x 



and thus 



\yy ~ "^ ,y 



v' 



V20 = ^ + 
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dLxx 


d(f) 


J^xx 


dz 


dz 
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dLxx 


d(p 


J-iXX 


dz 


dz 


d^<p 







dz^ 
2 dL, 



8.5.3 
8.5.4 
8.5.5 
8.5.6 

8.5.7 

8.5.8 



dz"^ Lxx dz dz 

It is easy to see that this leads to the correct Laplacian for flat space (i.e. put Lxx = f^^d 
and r = z where r is the usual radial coordinate and z is the proper distance measured along 
the radial axis). 



8.6. The ADM Constraints 

The Hamiltonian constraint 



= R + K^ -K^^K^,, (4.2.1^ 



can be evaluated directly on the lattice as each term on the right hand side is known at each 
vertex of the lattice. Furthermore, since the lattice is spherically symmetric the only terms 
which survive are those that contain Rxyxy, Rxzxz, Kxx and Kzz- This leads to 

U = -Tixyxy ~r ^ii'xzxz ~\~ -t^xx ' ^^xx^zz \^.Z.o) 

The standard form of the ADM momentum constraints are 

= i^l^-i^/|, (8.6.1) 
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where K = K^ ^. These equations require more care as they contain covariant derivatives. At 
the origin of a RNC frame we have g^p = diag(l, 1, 1) and = T^^, and thus the constraints 
may be reduced to 

= K,^ - K/^p (8.6.2) 

Each of the partial derivatives can be evaluated using the methods of the previous section 
(though modified for use on a two index tensor). The results are 

K,^ = K,y = (8.6.3) 

^^^dK^^dKyy^dK^ (8.6.4) 

dz dz dz 

1 dLxx , j^ j^ TZ \ '^^z 



K/,^ = -^ {2K,, - Kxx - Kyy) + -^ (8.6.6) 



One also finds that terms such as K^z^x 7^ even though Kxz = at the origin of each RNC 
cell. This fact has been used in the above results. 

The only non-trivial momentum equation is that ioT fi = z and this leads to 

= ^ (LxxKxx) _ j^^^^h^ (4 2.4) 

dz dz 

where we have used the fact that K^x = Kyy in our RNC frame. 

8.7. Bianchi Identities 

The Bianchi identities in a RNC frame are just 

= R^ual3,p + Rfivf3p,a + Rfivpa,f3 (8.7.1) 

In a spherically symmetric space there is only one non-trivial Bianchi identity, namely, 

U = -tixyxy,z ' J^xyyz,x + J^xyzx,y \P. I ■^) 

The only non-zero components of the Riemann tensor at the origin of a RNC cell are Rxzxz = 
Ryzyz and Rxyxy (and others obtained by standard symmetries in the indicies). However, 
like the calculations above for the momentum constraints, we find that many terms including 
Rxyyz,x and Rxyzx,y are not zero. This fact is a simple consequence of the mixing that occurs 
amongst the non-zero R^uap brought about by the rotation matrices. Rather than list all of 
the non-zero derivatives we shall list only those that we need for the above Bianchi identity. 
They are 



dRx 

dz 



T3 ^xyxy /Q y fjN 

J^xyxy,z — j_ yo.i.o) 
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1 dL 



J^xyyz,x — J , y^xzxz J^xyxy) [p. I A) 

^xyxz,y ^ T ~] {^xzxz J^xyxy) [o.t .o) 



which when substituted into the above equation leads to 



= ^^i%^^-i?....^ (3.2.4) 

az dz 



8.8. Non-uniform finite diff'erences 

By applying standard Taylor series expansions to a smooth function f{z) it is easy to derive 
the following second order accurate finite difference approximations 



U,^ Lj^^ -T -(^22: \ \ ^zz / \ ^ZZ 



d'^f 2 fr r ^ f f° I .g^g^2) 



dz iJ y^ ~\~ -Li ^^ \ Li.y^ ±j 



zz ' zz \ zz zz 



for a non-uniform lattice (where L^^ ^ L"^ are the lattice spacings). Centred finite differences 
are not appropriate for two simple reasons. First, we chose our initial data Lzz to be non- 
uniform. Second, even if we did choose an initially uniform lattice, the subsequent dynamics 
{dLzz/dt 7^ 0) would immediately produce a non-uniform lattice. 

These approximations are used at various places in the text (e.g. for dL^x/dz and d'^N/dz'^). 

The only exception to the above was in the discretisation of the Bianchi identities. This 
equation was approximated at the centre of the radial struts by a forward finite difference 
operator df /dz = (/""" — f°)/ L^^ and by setting Rxzxz = {Rxzxz + ^xzxz)/'^- 
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Figure captions 

Figure 1. An embedding diagram for a 2-dimensional slice of the spherically symmetric 
3-geometry of a Schwarzschild black hole. The lattice is constructed as a ladder, with one 
end on the throat and the other in the weak field region of the black hole. The two radial 
edges of the ladder are radial geodesies, while each rung is a geodesic segment of the full 
3- metric. Thus each rung is not confined to the 2-spheres (except at the throat). 

Figure 2. The ladder on which the numerical solution was built. Successive computational 
cells overlap by sharing two successive rungs of the ladder. In this diagram there are just three 
computational cells. The production runs employed 800 rungs and had an outer boundary 
set at approximately r = 200M. 

Figure 3. A typical computational cell. The coordinate frame has been oriented so that 
the z-axis points in the usual radial direction while the origin has been located so that the z 
coordinate of vertices (1) and (2) equals zero. The y— axis points directly into the page and 
thus it has been suppressed. 

Figure 4. The three neighbouring computational cells used in applying the Binachi iden- 
tities. For simplicity we have not drawn the middle rungs in each cell and we have drawn 
each geodesic segment as a straight line. Once again the y-axis has been suppressed. 

Figure 5—13. These figures show the evolution of the basic lattice data for t = 10m to 
t = 100m in steps of 10m and also from t = 100m to t = 1000m in steps of 100m. All of 
the figures display a smooth evolution with no signs of any instabilities. In each of these 
figures we can clearly see the stretching that occurs in the grid. On each curve we have 
used a diamond to mark the location of the apparent horizon. In many of the figures for 
< t < 100m we can clearly see the loss of resolution brought on by the grid stretching 
(e.g. the fall off in the plateau of K^x and in the decay in the sharp peak of K^z)- These 
effects are much more pronounced in the long term evolution, 100m < t < 1000m. Note 
that at these late times the radial legs near the apparent horizon have been stretched by 
almost a factor of 100 while the rungs have been shrunk be a factor of approximately 100. 
This is a severe change in shape of the grid and so its not surprising that the accuracy has 
been lost. Each of the plots in figures 5-13 were produced using the grid centred scheme on 
a Bernstein, Hobill and Smarr grid with 800 grid points. The plots for < t < 100m were 
restricted to a proper distance of 100m simply to better display the changes in the grid. At 
t = 100m the grid extends out to a proper distance of over 265m. 

Figure 14. This displays the curvature term Rxyxy at t = 100m for four different resolutions 
of 100, 200, 400 and 800 grid points. This clearly shows that the ability to maintain a 
fiat plateau behind the horizon is compromised when there is a loss of resolution near the 
apparent horizon. 

Figure 15. This shows the exponential collapse of the lapse at the throat. The work of Beig 
[24] shows that N{r = 0) ~ e thus a plot of InA^ versus t should be a straight line. For 
t = to t = 100m the line is very straight, while for longer times a slight bend does occur. 
The lapse at the throat at t = 100m is approximately 3.4 x 10^^^ while at t = 1000m is the 
lapse is of the order of 10~^^^. 



Figure 16. The size of tlie apparent liorizon from t = to t = 100m for four different 
resolutions of 100, 200, 400 and 800 grid points. 




Figure 1. The embedding diagram. 




Figure 2. The ladder lattice. 
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Figure 3. The typical computational cell. 
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Figure 4. Cells for the Bianchi identities. 



Figure 5. Leg lengths 
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Figure 6. Leg lengths 
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Figure 7. Riemann curvatures 
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Figure 8. Riemann curvatures 
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Figure 9. Extrinsic curvatures 
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Figure 10. Extrinsic curvatures 
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Figure 11. Hamiltonian constraint 
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Figure 12. Momentum constraint 
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Figure 13. Lapse function 
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Figure 14. Resolution 
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Figure 15. Collapse of the lapse 
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Figure 16. The apparent horizon 
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